function y = func_delta_from_p(p, m, CONSTS)
    
    k0    = CONSTS.k0;
    a     = CONSTS.a;
    eta   = CONSTS.eta;
    eps_a = CONSTS.eps_a;

    s       = func_s_from_p(p, CONSTS);
    q       = func_q_from_p(p, CONSTS);
    [n1 n2] = func_nk_from_p(p, CONSTS);
    K_m     = func_K_m_from_p(p, m, CONSTS);
    [J_m1 J_m2 J_m1_1 J_m2_1] = func_J_mk(q, m, CONSTS);
    
    S = k0*a*s;

    y = (n2-n1).*((p.^2)./eps_a).*((m^2)./((S).^4)) - (n2-n1).*(K_m).^2 + ...
        n1.*(eta/eps_a).*(J_m1_1).*(J_m2_1) - n2.*(eta/eps_a).*(J_m1).*(J_m2) - ...
        n2.*J_m1.*K_m + (eta/eps_a)*n1.*J_m1_1.*K_m + n1.*(J_m2_1).*(K_m) - ...
        (eta/eps_a)*n2.*(J_m2).*(K_m) + p.*(eta/eps_a).*(m./((S).^2)).* ...
        (-(J_m2_1) + (J_m1) + (J_m1_1) - (J_m2));

end

function y = func_K_m_from_p(p, m, CONSTS)

    k0 = CONSTS.k0;
    a = CONSTS.a;
    
    s = func_s_from_p(p, CONSTS);
    
    S = k0*a*s;

    y = (besselk(m+1, S, 0)./((S).*besselk(m, S, 0))) - ...
           (m./((S).^2));
       
end
